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Abstract. We present an "equation-free" multiscale approach to the simulation of unsteady 
diffusion in a random medium. The diffusivity of the medium is modeled as a random field with 
short correlation length, and the governing equations are cast in the form of stochastic differential 
equations. A detailed fine-scale computation of such a problem requires discretization and solution of 
a large system of equations, and can be prohibitively time-consuming. To circumvent this difficulty, 
we propose an equation-free approach, where the fine-scale computation is conducted only for a 
(small) fraction of the overall time. The evolution of a set of appropriately defined coarse-grained 
variables (observables) is evaluated during the fine-scale computation, and "projective integration" 
is used to accelerate the integration. The choice of these coarse variables is an important part of the 
approach: they are the coefficients of pointwise polynomial expansions of the random solutions. Such 
a choice of coarse variables allows us to reconstruct representative ensembles of fine-scale solutions 
with "correct" correlation structures, which is a key to algorithm efficiency. Numerical examples 
demonstrating accuracy and efficiency of the approach are presented. 

Key words, multiscale problem, diffusion in random media, stochastic modeling, equation-free. 

1. Introduction. This paper is devoted to numerical simulations of diffusion in 
a random medium whose material property, i.e., diffusivity (permeability, conductiv- 
ity), is characterized by small-scale, rough structures. This problem arises in the study 
of composite material properties, flow in porous media, etc (e.g. |51l31|). Direct, fully 
resolved computations of the governing equations in such media can be prohibitively 
time consuming, as the fine-scale structures require discretizations resulting to large 
degree-of-freedom calculations. Hence, there has been a growing interest in designing 
efficient alternative methods to solve the problems with the desired accuracy. 

The properties of such media are typically modeled as deterministic smooth func- 
tions superimposed with fast oscillatory components. One of the traditional ap- 
proaches is to derive the effective properties of such (heterogeneous) media - the 
so-called "homogenization" or "upscaling" techniques. These techniques are typically 
based on analytical asymptotic treatments and have been remarkably successful in 
several applications. Their applicability may, however, be restricted due to the as- 
sumptions that need to be made for the media (See, for example, 0EJ0). Numerical 
(as contrasted to analytical) approach to homogenization are typically based on build- 
ing multiscale basis functions into the spatial discretization. Methods along this line 
of approach can be found in ^3 1171 1181 131)1 and are the subject of active research. 

Alternatively, we can choose to model the medium properties as random fields, 
to account for our insufficient knowledge and/or measurement error. For example, 
field data indicate that the conductivity of many natural porous formations can be 
accurately described by a lognormal distribution (e.g. 9 ). The homogenization 
techniques developed for deterministic media can be generalized to random media, 
and a comprehensive review can be found in 

In addition to deriving equations describing an effective medium, as homogeniza- 
tion does, many efforts have been devoted to direct detailed simulations of random 
media. In this context, the corresponding governing equations, e.g. the diffusion 



* Department of Chemical Engineering, Princeton University, Princeton, NJ 08544. 
(xiu@dam . brown . edu) . 

' Department of Chemical Engineering, Princeton University, Princeton, NJ 08544. 
(yannisOprinceton. edu) . 



1 



2 



DONGBIN XIU AND IOANNIS KEVREKIDIS 



equation, Darcy's law, etc., are cast in the form of stochastic equations and solved 
as such directly. This approach further complicates the problem, since the governing 
equations are now defined in (much) higher dimensional spaces, composed of both 
the physical space and the space accounting for the parameterization of the medium 
randomness (the random space). The most straightforward numerical approach is 
the Monte Carlo method, (e.g. (HJ) where repetitive deterministic simulations are 
conducted for particular realizations of the random functions describing the medium 
properties, and what we are interested in here is the statistics of the solution. This 
approach, based on random sampling, can be computationally expensive because the 
convergence rate of the ensemble averages, e.g., mean solution, is relatively low (Monte 
Carlo simulations consisting of M realizations converge at a rate of 1/ vM-) There has 
been, therefore, a continuing interest in constructing non-sampling methods, which 
include perturbation methods [22], second-moment analysis [2E], stochastic Galerkin 
methods [691 131 ITS] . etc. Among them, the stochastic Galerkin methods, also called 
"generalized polynomial chaos" expansions, have been successful in many applica- 
tions, when the basis functions in the random space are appropriately chosen. In 
particular, when the solution is sufficiently smooth in the random space, stochastic 
Galerkin methods converge exponentially fast 001 CUD EH • 

However, for the problem we intend to study in this paper - random media 
with short correlation length - the stochastic Galerkin methods become inefficient. 
This is because the short correlation length will induce a higher-dimensional random 
space, subsequently larger number of equations to be solved - the so-called "curse- 
of-dimcnsionality" . We remark that such a difficulty exists for all the existing non- 
sampling methods. On the other hand, although Monte Carlo methods are (formally) 
immune from such a curse-of-dimensionality, their inherently slow convergence rate 
can hardly improve the overall efficiency. 

In this paper, we propose an equation-free, multiscale method to simulate dif- 
fusion in a random media with small-scale structure (short correlation length). The 
equation-free methods were first introduced in |88j and are designed to resolve multi- 
scale problems efficiently. Such methods solve the equations for the effective, coarse- 
grained behavior without obtaining them in closed form; the quantities required for 
these computations (residuals, action of Jacobians, time derivatives) are estimated by 
solving the microscopic/stochastic model with appropriately chosen initial conditions 
over short times (and, in certain cases, only parts of the spatial domain). The numer- 
ical results of these appropriately initialized short bursts of microscopic computations 
are used to estimate the rate-of-change (or other quantities of interest) of appro- 
priately defined observables: macroscopic variables characterizing the coarse-grained 
evolution. These rates are then used by projective integration to evolve the coarse- 
grained observables in time with (hopefully much) larger time steps (^] El OS])- 
Thus, the time-consuming microscopic solvers are used for only a (small) fraction of 
the overall time integration and no explicit knowledge of the equations governing the 
macroscopic variables is required (equation-free). This framework has been applied 
to a variety of problems, ranging from bifurcation analysis of complex systems to 
homogenization of periodic media [TT1 IHll2"Tl 1231 I2"gll2"§l 1551 I5B] . 

For the random media with rough structure considered in this paper, we em- 
ploy the Monte Carlo method as the "fine-scale" solver. An orthogonal polynomial 
expansion of the fine-scale solution is conducted (in principle) at every point in the 
physical space. Our coarse-grained observables are the first few expansion coefficients 
of such pointwise polynomial expansions on a relatively coarse grid; the key assump- 
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tion underlying our method is that in principle it should be possible to write a closed 
equation for these observables that can successfully describe the (long term) evolution 
of the solution statistics. 

The particular assumption (observation) in this paper is that, although each in- 
dividual realization of the solution is characterized by small spatial scales, induced by 
the small scales in the diffusivity field, ensemble solution averages are characterized by 
larger, "coarser" scales. That is, after possibly a short initial transient (relaxation), 
the ensemble solution averages are significantly smoother than individual realizations 
in space, and can therefore be accurately approximated with (hopefully significantly) 
fewer degrees of freedom (e.g. on a coarse mesh). This, in turn, implies that a closed 
equation exists (whether we can explicitly derive it or not) for these averages on a 
coarser mesh; this is precisely the equation that we will try to solve in this paper, 
without explicitly deriving it. 

The coefficients of the pointwise polynomial expansion of the random solutions 
are representative of such ensemble averages, and are observed to be smoother func- 
tions in space. We thus expect that they can indeed be represented on a coarse mesh. 
Since the explicit form of the governing equations for the evolution of such coefficients 
is unknown (to our best knowledge), we employ the equation- free framework to com- 
pute with it. In effect, we are trying to combine the simplicity of the Monte Carlo 
implementation with the strengths of (generalized) polynomial chaos representation: 
instead of deriving and discretizing the equations for the appropriate polynomial chaos 
coefficients via Galerkin expansion, we try to solve these equations through the design 
of "just enough" short computational experiments with the detailed direct solvers. To 
this end, the rate-of-change of these coarse-grained variables are estimated numeri- 
cally from short bursts of fine-scale computation, and propagated in time with larger 
steps via the projective integration technique. The advantage of the present definition 
of the coarse variables is that it allows us to reconstruct representative ensembles of 
fine-scale solutions with controlled accuracy. Numerical examples are presented to 
document the accuracy and efficiency of the new algorithm. 

The paper is organized as follows: in Section [2] we formulate the mathematical 
framework for diffusion in a random medium, and subsequently the multiscale prob- 
lem we will study. In Section |31 we present the details of the construction of our 
"equation-free" multiscale algorithm; in particular we focus on the "lifting" step": 
the construction of representative ensembles of fine scale solution realizations. Nu- 
merical examples are presented in Section^] and we conclude the paper with a general 
discussion in Section 

2. Unsteady diffusion equations in random media. In this section, we be- 
gin by presenting the mathematical framework for diffusion in a random medium. We 
then formulate the multiscale problem that we will study and discuss the difficulties 
in solving it efficiently. 

2.1. Formulations for random diffusion. Let D 6 M d , d = I, 2, 3 be a bounded 
polygonal domain, J = (0, T] 6 M. + = (0, oo) for some fixed time T > 0, and (f2, J 7 , P) 
be a complete probability space. Here fl is the set of outcomes, T C 2 is the 
er— algebra of events and P : T — ► [0, 1] is a probability measure. Let Dt = D x J, 
and we study the following random diffusion equation: find a stochastic function, 
u:Ox Dt —* R, such that for P— almost everywhere u> G CI, the following equation 
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holds: 

u t (oj,-) ~W ■ (k(u,x)Wu(uj,-)) = f(uj,-), in D T) 

u(uj,-) = 0, on dD X [0,T], (2.1) 

u(lu,-) = uq(lo,x), on D x {t = 0}, 

where u is the unknown and Ut — du/dt its time derivative. k,uq : f2 x D — > 
R, and / : f2 x £> T — > R are the known stochastic functions with continuous and 
bounded covariance functions. Denote by B(A) the Borel a— algebra generated by the 
open subsets of A, then k and uq are assumed to be measurable with the a— algebra 
(T ® B(D)) and / with (T ® B(D T )). 

The following assumptions are made on the input stochastic data: 

1. ft is bounded and uniformly coercive, i.e. 

3 K min , ft max e (0, +oo) : P (w G CI : k(u>,x) G [ftmin, K max ],Vx € D) = 1. 

(2.2) 

Also, ft has a uniformly bounded and continuous first derivative, i.e. there 
exists a real deterministic constant C such that 

P(we!l: ft(w, •) G C 1 ^) and max s |V x k(w, -)| < C) = 1. (2.3) 

2. f e L 2 (n)(g> L 2 (D T ), i.e. 

/ / / f 2 (iu,x,t)dxdtdP{uj) < +oo. (2.4) 
ij Jd 

3. w G L 2 (n)®L 2 (D), i.e. 

UQ(u>,x)dxdP(ui) < +oo. (2-5) 

2.2. Finite-dimensional noise and variational form. For the problem (|2.1Jl 
to be practically solvable numerically, it should be possible to reduce the infinite- 
dimensional probability space to a finite-dimensional space. This can be accomplished 
by characterizing the probability space by a finite number of random variables. Such a 
procedure, termed as the "finite-dimensional noise assumption" in [3], is often achieved 
via a certain type of decomposition which can approximate the target random process 
with desired accuracy. One of the choices is the Karhunen-Loeve type expansion 
|27|. which is based on the spectral decomposition of the covariance function of the 
input random process, (e.g. |15l I39| ). Following a decomposition and assuming that 
the random inputs can be characterized by N random variables, we can rewrite the 
random inputs in the following abstract form, 

k(u,x) = k(Yi(ui), . . .,Y N (ui),x), and f(u>,x,t) = f(Yi(ui), . . . ,Y N (w),x,t), (2.6) 

where N > 1 is a finite integer and {Y n }^ =1 are real random variables with zero 
mean value, unit variance, and their images T n ^ = Y n (Q,) are bounded intervals 
in R for n = 1, . . . , N. Moreover, we assume that each Y n has a density function 
Pn ■ ^n,N — * R + for n — 1,...,N, and denote p(y),Vy G T the joint probability 
density of (Yi, . . . , Y n ) and V = Jl^Li ^ n ,N C R w the support of such density. The 
expectation operator is subsequently defined as E(/) = L f(y)p(y)dy. Note when 
random variables Y n , n = 1 . . . , N are independent, we have p(y) = n«=i Pn(yn),Vy G 

r. 
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After the finite-dimensional characterization of the random inputs l|2.6[l , the un- 
steady diffusion equation l|2.1[l can be expressed in the following strong form, 

u t (y,x,t) - V- (n(y,x)Vu(y,x,t)) = f(y,x,t), V(y,x,t) G r x D x J, 

u(y,x,t) = 0, y(y,x,t) G T x dD x [0,T], 

u{y,x,0) = u (y,x), V(y,x)ETxD. 

(2-7) 

Often we seek its weak solution satisfying the following variational form: find 
u e L 2 p (T) <x> L 2 (0, T; H^(D)) with u t G L 2 p (T) <g) L 2 (0, T; H-\D)) such that 



where 



and 



!„(«*, =I P (f,v), Vv € L 2 p (F) ® H%{D), 

u(t = 0) = Uq 



l p (v,w)= / p(y) / v(y,x)w(y,x)dxdy 7 
Jr Jd 



JC p (v,w;k)= / p(y) / n(y,x)Vv(y,x) ■ Vw(y,x)dxdy. 



(2.8) 



Note that problem 1)2. 7fl or (|2.8I) becomes a (AT + d) dimensional problem, where 
d is the dimensionality of the physical space D and is the dimensionality of the 
random space T. 

2.3. Formulations for multiscale problems. In this section we formulate 
the multiscale problem associated with the stochastic diffusion equation l|2.1|l . In 
particular, we consider the problem where the random inputs, e.g. diffusivity k, 
have very small correlation length 1 N C 1, compared to the (macroscopic) domain 
of interest Dt- For notational convenience, hereafter we restrict our exposition to 
problems with only k being the random input, and study, for P-almost everywhere 
uj G ft, 

^(u;,.)=V- [«(o;,J)w(uv)] +f(x), mD T , (2.9) 



u £ (uj,-) = 0, ondDx[0,T], (2.10) 



u e (uj,-) =u e (x) on D x {t = 0}. (2.11) 

Here we have assumed that the diffusivity k satisfies the conditions (|2.2|l and i|2.3|) . 
and is a homogeneous random field with a short correlation length l K ~ 0(e) <§; O(l), 



C K (x, y) = C(\x - y\/l K ), l K ~ 0(e) « 0(1), (x, y) G D x D, (2.12) 

where C v (x,y) = E [(v(ui, x) — E[v(ui, x)])(v(ui, y) — E[v(cj, y)])] is the two-point co- 
variance function and l K is the correlation length. Again we characterize the diffusiv- 
ity field K(u},x/e) by N independent random variables as in (|2.6|) . Hence, problem 
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(|2.9|) is in (TV + d) dimensions, and we can formulate it in the strong and weak forms 
similar to (|2.7[1 and (|2.8|) . respectively. 

The discretization in the spatial domain D C M. d can be conducted via any stan- 
dard technique, e.g., finite difference, finite elements, etc., with a maximum mesh 
spacing parameter S > 0. To fully resolve i|2.9|) . we need to employ a fine mesh 
with 6 < e. From a numerical point of view, very small mesh spacing S often results 
in restrictions on the size of time steps of numerical schemes, and such restrictions 
are particularly severe for explicit schemes. Hence, a fine-scale computation of (|2.9[1 
requires computations with very small time steps on a very fine mesh. 

The discretizations in the TV— dimensional random space T can be conducted in 
different ways. The recently developed stochastic Galerkin methods, or generalized 
polynomial chaos, are extensions of the classical polynomial chaos which is based 
on the Wiener-Hermite expansion f5 . These extensions include the non-Hcrmite 
global orthogonal polynomial expansion from the Askey scheme 001113021) piecewise 
polynomial expansions El EJ , and wavelet basis [23 E] • The stochastic Galerkin 
methods have fast convergence as the polynomial order is increased. In fact, un- 
der sufficient regularity requirements, exponential convergence has been proved for 
stochastic elliptic equations in [3] and shown numerically for various stochastic equa- 
tions in HOI EHlllIl- The total number of expansion terms, K , however, depends not 
only on the order of polynomials but also the dimensionality TV of the random space. 
When TV 3> 1 is very large, K increases fast with increasing order of polynomials. 
This significantly reduces the convergence rate with respect to the number of expan- 
sion terms for stochastic Galerkin methods. To this end, it may be necessary to resort 
to the Monte Carlo method, as its convergence rate, 1 / \f~M where M is the number 
of realizations, albeit slower, is independent of the value of TV. 

The number of random variables, TV, used to represent the random process k(lu, x) 
depends on, among other factors, the correlation length of k. Although one may 
choose different decomposition methods, in general, the value of TV is inversely pro- 
portional to the correlation length l K . For the problem we consider here, l K <C 0(1) 
implies TV 3> 1, and problem (|2.9|) is in a high-dimensional random space. Subse- 
quently, the fine-scale computation of (|2.9I) requires a large number of discretization 
terms in the TV-dimensional random space T (by either a large number of expansion 
terms from a stochastic Galerkin method at a given order, or a large number of real- 
izations from a Monte Carlo method) , a fine mesh in the physical space D to resolve 
the small scales induced by l K <C O(l), and very small time steps in the time domain 
J. Such computations can be prohibitively time-consuming. 

3. An equation-free multiscale method. In this section we present an equation- 
free multiscale method for the integration of equation (|2.9|l ; other tasks (such as fixed 
point algorithms for its stationary states) can also be formulated in an equation-free 
context (see the discussion in Section EJ). The key feature of the method is that the 
costly fine-scale computations of (|2.9|l are only conducted for a small fraction of the 
total integration time. During the fine-scale computations, the rate of change of ap- 
propriate coarse-scale variables is estimated numerically and subsequently represented 
on a coarser mesh and integrated in time with large time steps. The choice of such 
coarse variables is based upon the following assumption (observation) : although each 
individual realization of the solution of H2.9[) is characterized by small spatial scales, 
the ensemble averages, e.g. moments, of the solutions are significantly smoother in 
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space (characterized by much larger scales), i.e. 

= U g (x,t), VffGC, (3.1) 

where C denotes a set of smooth functions. Variables U g (x, t) are the "coarse-grained" 
variables, and we will describe in detail their construction in the following Section. 
Hereafter, we will drop the superscript e of the fine-scale variables u 6 . 
Equation l|2.9[l defines an evolutionary process 

— (w,x,t)=r(u), (w,x)eQxD, (3.2) 

characterized by a solution operator {s(t)}, which forms a semigroup u(-, t) — s(t)u{-, 0) 
in t £ J. We will assume that the set of properly defined coarse-scale variables from 
(13.11) satisfy, possibly after a short transient (relaxation) period, closed differential 
equations 

BTJ 

-^-(x,t) = R(U g ), (x,t)£D T . (3.3) 

Note that typically there are multiple number of coarse variables. Subsequently, U g 
is a vector field and 1)3. 3|1 is a system of equations. We also remark that the explicit 
knowledge of the equation 1)3.3(1 may not be easy to obtain, or may be too complicated 
to be of any practical use if it were known. 

The general procedure for the equation-free projective integration methods (see, 
for example, ^JI^EOI) over one global time step At, starting at t = t™ and ending 
at t = t n+1 , consists of the following key components: 

• a "restriction" operator V to evaluate the coarse-grained variables from the 
ensemble of fine-scale computations, i.e. U g — Vu, and a "lifting" operator 
Q to construct the representative ensemble of fine-scale solutions from the 
coarse-scale variables, i.e. u = QU g ; 

• rif > 1 steps of fine-scale computations of (|3.2(l with a small time step 5t, 
where we will define At/ = rifSt and an intermediate time level t™ = t n + Atf; 

• one step of coarse projective integration of the coarse-scale equation (|3.3|l 
with a time step of the size At c = n c St, n c > 1. 

Since At c is associated with the time scale of the coarse variables U g defined in 1(3. 1() , 
we usually have n c » 1. The global time step is At = t n+1 — t™ = At/ + At c = 
(rif + n c )8t. Figure l3~D is a graphical illustration of the notations. 

Specifically, for the multiscalc diffusion problem in a random medium described 
by 1)2.9(1 . the equation- free integration over one global time step At consists of the 
following steps: Given a fine and coarse computational mesh, and U™(x) = U g (x,t n ) 
on the coarse mesh, 

1. Lifting (or Reconstruction): Generate an ensemble of random solutions u n (uj, x) = 
u(uj,x,t n ) on the fine mesh. 

2. Fine-scale computation: Fully resolve l|2. 9J) by using the fine mesh in D, a 
small time step St in J, and an appropriate method in T (e.g., the Monte 
Carlo method). Such a fine-scale integration is conducted only for a short 
period of time, from t™ to the intermediate time t™, i.e. u(-,t) — s(t)u(t n ), 
for t™ < t < t™ = t n + nfSt. Here rif > 1 such that nfStf ~ tn ^ tjvf, where 
tn is the local relaxation time of the fine-scale process and it is assumed to 
be much shorter than tM, the typical time scales of the coarse-scale process 

63> . 
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Projective integration 

Exact 




Fig. 3.1. Sketch of the multiscale equation-free integration over one global time step. 



3. Restriction: Evaluate the coarse variables U g {t) defined in (|3.1fl on the coarse 
mesh, for t n <t< t n c . 

4. Coarse-scale integration: Estimate the time derivatives of the coarse variables 
U g at t = i™ and integrate the coarse-scale equations (|3.3|) to t n+1 via the 
projective integration method on the coarse mesh, with a time step At c — 
n c St = t n+1 - t n c . Here At c ~ t M > t R . 

We now present a detailed description of each of the steps, starting with the more 
straightforward step - the fine-scale computation. 

3.1. Fine-scale computation. The objective of the fine-scale computation is 
to fully resolve (|2.9(l . To this end, any conventional spatial and temporal discretization 
scheme can be employed, e.g. finite difference or finite elements. Since the dimen- 
sionality is high ((N + d)-dimcnsional), we employ a Monte Carlo simulation (MCS) 
in the random space T. Here we illustrate the formulation of a Monte Carlo finite 
element method (MCFEM). 

Denote Xf C Hq(D), flcl^a family of piecewise linear finite clement approxi- 
mation spaces, with a maximum mesh spacing parameter 5 > 0. This is the fine mesh, 
as we choose 5 < 0(e) <C O(l) to fully resolve all spatial scales. We shall assume all 
the standard assumptions on the finite element triangulation, and its approximation 
estimate, i.e. 

min \\v-x\\<CS\\v\\ HHD) , Vi> G H 2 (D) n Hq(D), (3.4) 

where C > is a constant independent of v and 8. 
The MCFEM formulation for is as following: 

• Prescribe the number of realizations, M, and a piecewise linear finite element 
space on D, Xf, defined as above. 

• For each j = 1, . . . , M, sample i.i.d. realizations of the diffusivity , •) and 
find the corresponding approximation u$(ujj,-) G L 2 (0,T; Xf) with ^jf- £ 
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L 2 (0,T; {Xf)') such that 

r^( WjV ),x) + J K(u> j ,-)Vu s (u j ,-)-v x dx = (f, x )s, v x exl,teJ, 

(3.5) 

where (■, ■)$ is the usual inner product in Xf. 
• Process the solution ensemble to generate statistics, e.g. E(u) — J2jLi u s {^j > ' 
For more detailed discussion on the stochastic finite element spaces for elliptic 
problems, see for numerical examples and implementations of stochastic Galerkin 
methods for steady/unsteady diffusion equations, see |391 142j . 

3.2. Restriction and Lifting. The restriction from the fine-scale variables u 
to coarse-scale variables U g consists of two steps: "random restriction" and "spatial 
restriction". First, the fine-scale variables u are averaged to U g in the random space 
according to l|3.1J) (random restriction); then the coarse variables U g are further re- 
stricted from the fine mesh to the coarse mesh (spatial restriction) justified by the 
assumption/observation that they are smoother. The lifting procedure is the reverse 
of the restriction. We now describe the details of the two procedures in both the 
random space Y and the physical space D. 

3.2.1. Operations in random space. The solution of the fine-scale computa- 
tion via MCFEM in section mi or other effective methods, generates an ensemble of 
M realizations of the random solution ug at any x € Xf and t E J. For any fixed 
(x,t), u$(uj, ■) is a random variable, and we seek to represent such a random variable 
by an orthogonal polynomial approximation and define the expansion coefficients as 
our coarse-grained observables (variables). Hereafter we drop the subscript S in us, 
the fine-scale numerical solution of u, and seek to approximate it by I u u for any fixed 
(x,t), 

K 

u(ui, x, t) ~ 1 u u(lu, x, 

k=0 

where {$fc(£M)HL o is a set of orthogonal polynomials {$} in term of random vari- 
able The expansion coefficients are determined by 



U k = / n u(uv)**(£(w))dPH 



E <I>; 



T E [ttfo •)**(£("))], k = 0,---,K, ( 37 ) 



where the orthogonality of the basis functions has been used. Such a representation of 
random variables is commonly used in practice (cf. 15, 40 39 ). The correspondence 
between the type of orthogonal polynomials {$} and the type of random variable £(w) 
includes Hermite-Gaussian, Jacobi-beta, Laguerre-gamma, etc (see @U] for details). 
The convergence of such orthogonal polynomial expansions is assumed to be of the 
form 

\\u(w,')-X u u(w r )\\v {a) =O(K-i'')^0, K^oo, 7w >0, (3.8) 

where we have assumed the convergence rate scales as K^ 1 ^ for some positive number 
7^ > 0, which depends on the smoothness of u(cu). We remark that a complete 
theoretical analysis on the convergence of different basis remains an open issue. For 
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numerical examples of the approximations of a random variable via different sets of 
basis, see |IU] . 

The expansion coefficients {U} k=1 are the ensemble averages of it, as defined 
in H3.7jl . and under assumption (|3.1[1 . they become the coarse variables with larger 
spatial scale, i.e., with smoother profiles in the physical space D. For instance, the 
coefficient Uo(x,t) is the mean field of u and is in general smooth. 

The finite-term polynomial approximation (|3.6|) defines two operations between 
u and {U k }, i.e., the "restriction" operator in random space Vu, such that, 

{#*(•)}£=<> =n.u(<«v)> (3-9) 
and the "lifting" operator such that 

l u u(u J ,-) = Q u) {U k (-)}^ , (3.10) 

where operation T w is accomplished by (|3.7|) , and Q u by generating random samples 
of in equation (|3.6|) . Obviously both V u and Q w are linear operators, QJP^ = I u 
and T^T^ = T^. 

We remark that the expansion Ij3.6|l is different from the traditional polynomial 
chaos expansion. Expansion l|3.6|l is a pointwise approximation at fixed locations in 
(x,t), and hence only requires a one-dimensional (in the random space) polynomial 
basis of {$fc(£(w))} where the random variable £(oj) associated with the basis is dif- 
ferent at different locations in the physical space. On the other hand, the traditional 
polynomial chaos expansion is written in the full N— dimensional random space for 
all locations of {x,t). While operators (|3.9|1 and l|3.10() are one-dimensional in the 
random space, they do not offer us an easy way to obtain the governing equations for 
the coarse variables U, shown as in i|3~3|) . (On the other hand, we can readily derive 
the governing equations in the N— dimensional random space via a Galerkin method 
if the traditional polynomial chaos expansion is employed.) We will show in the next 
section that we can circumvent the difficulty of not having the governing equations 
by using the "equation-free" approach. 

3.2.2. Operations in physical space. Since we have assumed the coarse vari- 
ables {U k (x,t)} are smooth in space, they can be accurately represented on a coarse 
mesh X^ C D, whose maximum mesh spacing A>i. Such a representation can be 
expressed as, e.g., in terms of polynomial approximations in the physical space D, 

L 

U k (x,t) ~I x U k (x,t) =J2u k!l (t)<t> k (x), (3.11) 
i=i 

where {(^(x)}^ are the basis functions in X^, and the expansion coefficients are 

defined as U k> i = (U k , 4>i)A/{4>h <^z)a> Vfc. Here (•, -)a is the usual inner product in 
X^ and we have assumed, for notational convenience, that the bases are orthogonal. 
The completeness of such bases yields 

\\U k (x,-)-X x U k (x,-)\\ xi =0{L-^-*) ^0, L^co, lx . k >0,Vfc, (3.12) 

where || • \\x^ ^ s an appropriate norm in X^ and ^ x ,k > quantifies the convergence 
rate, which depends on the smoothness of the underlying function U k . 

Similarly, we can define two operators between U and {£/;}: the restriction oper- 
ator V x in the physical space D such that 



{Uk,l(')}='P*{Uk(x,-)}, k = 0,--- ,K and I = !,-■■ ,L, (3.13) 
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and the lifting operator Q x such that 

T x {U k (x,-)} = Q x {U k ,i(-)}, fe = 0, and / = I,- •■ ,i. (3.14) 

Clearly, we have Q, X V X — T-x and X X X X = X x . 

3.2.3. Global restriction and lifting. The global restriction operator V and 
lifting operator Q are thus defined as 

V = V X V U such that {U K i{t)} =Vu{uj,x,t)- (3.15) 

and 

Q = QujQx such that Xu(cj,x,t) = Q{U k ,i(t)}, (3.16) 

where the global approximation operator X is defined as 

X = QV = Q^Q x V x Vu = QJLxPu. (3.17) 

Note the above operators are defined for to = {u>j}jL 1 G ft, k — {0, ■ • • , K}, I — 
{1, • • • , L}. A remarkable property of the operator X is that 

Ilii-Xid] -> 0, K,L ^oo, (3.18) 

where the norm || • || is defined in the tensor product space of L 2 (fl) and the appro- 
priate space of D that defines the norm in l|3.12|l . Such a property ensures that, by 
restricting from the fine-scale solution ensemble to the coarse-scale variables (operator 
■p) and lifting back to the fine-scale (operator Q), we can reconstruct a representative 
ensemble of fine-scale solutions with controllable accuracy. 

3.2.4. Implementation of Restriction and Lifting. In numerical simula- 
tions, the ensemble of solutions obtained by the Monte Carlo method in Section l3~Tl 
{u(-, LJj^jL-y, are first restricted in random space by Ij3.9|l : the resulting variables 
{Uk{x)} on the fine mesh is further restricted to the coarse mesh by V x (|3.13() . 

The random restriction is accomplished by l|3.7|l , where the integral is written 
in a formal way as u and £ in the integrand do not in general have the same probability 
measure. To integrate (|3.7f) . we transform both the random variables u and £ to 
a uniform variable 9 £ J7(0, 1) via their cumulative density function (CDF), i.e., 
6 = F(u) = G(£), where F and G and the CDF of u and £, respectively. Hence, 

u = F-\9), ^G-^O), (3.19) 

and 1)3. 7|l can be written as 

1 f , 1 - 1 



Uk = f , 2 , p / u(w,.)* k (Z(w))dP(u) = / F-\e)3> k {G- l {6))d6. 

(3.20) 

The resulting integral in the bounded domain (0, 1) can be readily integrated using 
quadrature rules with sufficient accuracy. For more details on such integrations and 
numerical examples, see |4()j . 

During the lifting procedure H3.16J) . the coarse variables on the coarse mesh are 
first "lifted" to the fine mesh via Q x (I3.14|) (effectively, interpolated); the random 
lifting operation Q u is then conducted to generate an ensemble of fine-scale solutions 



12 



DONGBIN XIU AND IOANNIS KEVREKIDIS 



on the fine mesh. The spatial lifting Q x is accomplished through l|3.11|l . and the 
random lifting Q w by (|3.6|l . In (|3.(j|) . an ensemble of realizations of the random 
variable {^(wj)}^^ are generated to, in turn, generate {u^j)}^^. The {^(wj)}fLi 
is generated via (|3.19|) by using the same set of uniform random variable {6{uij)}Y=i 
resulted from the CDF of {u^j)}^^, i.e., 6 — F(u). By using the same set of 9, the 
lifted solution {Xu^j)}^^ will have the same correlation structure as {u^)}*^. 
This is a key step for efficient numerical computations, as it prescribes the "right" 
correlations between a particular realization of the medium and the corresponding 
solution for this medium. 

3.3. Coarse-scale integration. As pointed out in section 15.2. II although the 
pointwise expansion (|3.6I) only utilizes one-dimensional expansions in random space at 
fixed locations in physical space, it is unclear, to the authors' best knowledge, how the 
governing equations for the coarse-variables {Uk,i} should be, i.e. the right-hand-side 
of 

-^L = R ktl (U), Vfce [0,K], le [1,L\. (3.21) 

is unknown. To circumvent the difficulty, we employ the "equation-free" approach 
where explicit knowledge of these governing equations is not needed. This procedure 
is as following: 

• Evaluate the coarse variables {Uk,i(t)} = Vu(-,t) from the fine-scale compu- 
tation, for t" <t<t n c =t n +nf5t. 

• Approximate the RHS of Ij3.21jl at t = i.e. 

R k A(t n c ) = f^a 3 U k At 3 ) = ^(t:) + 0(St J f), ke[0,K], le[l,L], 

3=0 

(3.22) 

where 1 < n e < rif, tj = — jSt, and J/ denotes the order of the approx- 
imation, {ay}™^ IS a se t °f consistent coefficients such that J2 a j v (tj) = 
dv/dt{t^) + 0{St J f). 

• Once the RHS of l|3.21|l is estimated numerically, 13.21fl is integrated forward 
in time for one step on a larger time step. For example, given coarse time step 
of size At c = n c -St with n c > 1, such that t n+1 = t n + At c = t n + (nf + n c )8t, 
the Euler forward integrator takes the form 

Ukt 1 = UkAtc) + &tc-R k Atc) + o(At 2 c ), fee [o,K], le (3.23) 

The estimation of derivatives (|3.22|) may suffer from numerical oscillations, especially 
for systems that exhibit noisy behavior at microscopic scales, e.g., molecular dynamics. 
In this case, certain smoothing techniques such as a least-sqaure fit may be used to 
alleviate the problem ^1]. (Such is not the case in this paper.) The integration of 
coarse variables l|3.23|l is a simple Euler forward scheme, which can lead to relatively 
large errors when At c is large. Other integration schemes, such as higher-order single 
step methods or multi-step methods, can be used for their improved accuracy and/or 
better stability properties [111 I12| . A complete analysis (stability, accuracy, etc.) 
of the equation-free method to stochastic equations is still lacking, and is beyond 
the scope of the current paper. For an error analysis of equation-free method to a 
deterministic system (flow simulation), see |37|. 
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3.4. Computational Complexity. Let us denote by Nf ~ 0{8 d ) the number 
of degrees-of-freedom (DOF) of the fine mesh Xf, and iV c - 0(A~ d ) the DOF of the 
coarse mesh X^. During the coarse integration step, the fine-scale computation by 
M realizations of Monte Carlo simulations is effectively reduced to a problem of the 
evolution of (K + 1) local polynomial expansion coefficients l|3.6(l on the coarse mesh 
obtained by the global restriction operator V = VxVw The reduction in computa- 
tional complexity can be illustrated as 

R MxN f J^^K+lJxNf JP^ R (K+l)xJV c _ ( 3 24 ) 

Thus, to march problem (|2.9() over a global time step At — (nj + n c )5t, the par- 
ticular equation- free projective integration algorithm, which consists of rif steps of 
fine-scale computations and one-step of coarse integration, needs to solve a problem 
of complexity 

C c ~ n f x R M * N f + r(A'+i)xa^ (3 25) 

On the other hand, the complexity of the full-scale MCFEM over the same time 
interval At is 

C f ~ (n/ + n c ) xR MxN f (3.26) 

For the multiscale problem considered in this paper, we have Nf ^> N c , 0(1) ~ K <C 
M, and n c 3> nf. Thus, roughly speaking, C//C c ~ (l+n c /n/) 3> 1. We remark that 
such an estimate is rather crude and the actual computational efficency is problem 
dependent. 

4. Numerical Results. In this section we present numerical results on i|2.9|) in 

one spatial dimension x £ [0,1]. The restriction and lifting operations in physical 
spaces, as shown in Section 13.2.21 come from standard approximation theory, and 
here we focus on the properties of the random restriction and the random lifting. We 
remark that the method extends trivially to multi-dimensional physical spaces. We 
assume that k(u>, x) is a Gaussian random field with unit mean value, i.e., Eu(w, x) = 
1, and employ the Hermite polynomials in random space to represent the random field 
in 13.6fl . We employ the Gaussian random field model because of its ease of numerical 
generation. (Generation of non-Gaussian random fields is still an active research 
area.) From a mathematical point of view, Gaussian models are inappropriate for 
the diffusivity fields as they allow negative values with non-zero probability, and thus 
violate the uniform coercivity assumption l|2.2|) . In practice, however, such negative 
values are rare, especially when the variance is small, and we can neglect such negative 
values if they occur in the random realizations. (In the numerical examples below, 
negative values never appear.) 

We assume that the Gaussian random field has an exponential covariance function, 
i.e. C K {x,y) — exp(|a; — y\)/l K . Such a correlation function can be generated from 
a first-order Markov process, and has been used extensively in the literature. All 
MCFEM computations are conducted by a linear finite element method with Euler 
forward integration, and M = 1,000 realizations are used. The forcing term in (|2.9|l 
is fixed at f{x) = —2, and zero Dirichlet boundary conditions are imposed at x = 
and i = l. The restriction in physical space l|3.11|l is conducted on a set of Jacobi 
polynomial basis (see \HJ\, Chapter 2). 
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4.1. Accuracy. The complete error analysis of the present multiscale method 
remains an open issue. In this section, we conduct numerical experiments to document 
the error convergence. In the first example, we set l K = 0.1 and use 40 linear elements 
(S = 0.025). Fifth-order Hermite expansion [K = 5) is used for the random restriction 
(|3.6(l . and fifth-order Jacobi basis (L = 5) for the spatial restriction (|3.11|l . The time 
step for the fine-scale computation is St = 0.001 and it is conducted for n/ — 20 
steps. The coarse integration has a time step At c = 0.08 (i.e. n c — 80), so that 
the global time step is At = (rif + n c )Atf =0.1. For this illustrative example, the 
computational speed-up is modest (4 ~ 5 times), because the separation of scales is 
modest. 

Figure FTTI shows the stochastic solution profile (mean and standard deviation) at 
time T = 1. Good agreements are obtained between the full-scale MC simulation and 
the equation-free multiscale method. 




Fig. 4.1. Solution profile at T = 1. Left: mean solution, Right: standard deviation. 

To study the error contributions from different factors, we conducted a series 
of computations with varying parameters. In Figure 14.21 the L°° errors in mean 
and standard deviation (STD) are shown. These computations have fixed values of 
Atf = 0.05, K = 3, L = 4 and varying time steps At c of the coarse integration. We 
observe that the errors decrease as the size of time steps for the coarse integration 
decreases. 

To examine the error convergence with respect to the orders of approximation in 
random space (parameter K) and physical space (parameter L), we fix the time steps 
of integrations Atf = At c = 0.05. The size of coarse integration At c is sufficiently 
small such that the the temporal errors are subdominant (e.g., O(10~ 6 ) for the mean 
as shown from Fig. 14.2(1 . In Fig. 14.31 the errors with increasing order of the Hermite 
approximations (K) in the random restriction are shown. Fourth-order (L = 4) 
Jacobi basis is used in the physical space, so that the errors from spatial restriction 
V x are subdominant. It can be seen that as the order K of the random restriction 
increases, the errors in standard deviation decrease as expected. The mean solution is 
well-resolved by even the first-order Hermite expansion (K — 1) and its errors remains 
at the O(10~ 6 ) level. 

We then fix the order of the Hermite approximation in the random space at third- 
order, i.e., K — 3, and vary L. Again, At c = Atf — 0.05 is sufficiently small. In Fig. 
14.41 it can be seen that the error in the mean solution quickly reaches a saturation 
level of O(10 -6 ) at second-order L = 2. This is consistent with the result in Fig. 14.31 
The error in STD keeps decreasing with increasing order of L. 
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Fig. 4.2. Errors in mean and standard deviation (STD) versus the step size of the coarse 
integration. 



-0- mean 
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Fig. 4.3. Errors in mean and standard deviation (STD) versus the order of Hermite approxi- 
mation of the coarse variables. 



As discussed in Section I3.2.3I the current implementations of the restriction op- 
erator V and the lifting operator Q allow us to reconstruct representative ensembles 
of fine-scale solutions by restricting them to the coarse scale first and then lifting 
back to the fine scale, i.e. X = QV is an approximation operator. To examine 
the properties of the global approximation operator X l|3.17[) . we plot Au(uj,x,t) = 
u(u>, x, t) — Xu(oj, x, t) at an arbitrary chosen time t — 0.1. On the left of Figure ETBI 
several realizations of such errors (randomly chosen from the M = 1, 000 realizations) 
are shown. It can be seen that the errors are bounded within a small range of the 
same order of the spatial error (O(10 -3 )). The CDF of the random solution u at 
the center of the domain [x = 0.5) is shown on the right of Fig. 14.51 Again we see 
excellent agreement between the probability distribution of the target u and that of 
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Fig. 4.4. Errors in mean and standard deviation (STD) versus the order of spatial representa- 
tion of the coarse variables. 



its lifting Xu. 




To further examine the error of I, we plot on the left of Fig. 14.61 the point-wise 
error of Au(u),x,t) at x = 0.5, t = 0.1 for all realizations u>j,j £ [1, 1,000]. We ob- 
serve that, except at a few discrete points, which belong to a set with arguably zero 
measure, the errors are bounded in a very small interval of order O(10 -3 ). On the 
right of Fig. 14.61 we plot the path-wise correspondence of u(u>j) vs. lu(uij) for all 
j = 1, • • • ,1, 000. It is seen that the data collapse on y = x where the exact corre- 
spondence should be. These results confirm that our X is indeed an approximation 
operator which allows us to reconstruct the fine-scale solution ensemble with built in 
desired correlation structures from the computed coarse solutions. To achieve this, it 
is important that one uses, during the lifting step (|3.6(l . the same uniform random 
variable values determined by the numerical solutions at the restriction stage as de- 
scribed in Section ?6. 2. 41 If, however, a random variable £ is used without maintaining 
the correct correlation structure between the medium realization and the solution in 
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this medium, the lifted solutions will not be properly correlated to the true solutions, 
even if the ensemble is constructed to maintain the same distribution. Figure 14.71 
shows such an example. Again, this is the numerical solution at t = 0.1 and we plot 
the realizations at x = 0.5. Here the numerical solutions are lifted by using arbitrarily 
generated Gaussian random variables £ in (|3.6(l . We observe that although the solu- 
tion has the same distribution as the MCS solutions (Fig. 14.71 left), it is completely 
uncorrelated to the true solution as shown on the right of Fig. 14.71 The ability to 
maintain good correlation structure in the lifting procedure is a distinctive feature of 
our method. This is different from the conventional lifting procedures, whose recon- 
structed solutions are rather arbitrary (to a certain degree) and a certain constraining 
procedure or a relaxation integration is needed to "heal" the lifted solution ensemble 

HIT 3 




Fig. 4.6. Comparison of random solution (u) and its lifting (lu) at t = 0.1, x = 0.5. Left: 
point-wise error of \u(u>j) — Iu(u)j)\ for all = !,...,!, 000. Right: u(u>j) vs. Ju(wj). 
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Fig. 4.7. Comparison of random solution (u) and its lifting (lu) without correct correlation 
structure at t = 0.1, x = 0.5. Left: CDF of u and lu; Right: u(ujj) vs. Xu(uij) for j = 1,. . . , 1, 000. 

4.2. Efficiency. In the second example, we prescribe a diffusivity field with 
relatively small correlation length l K = 0.01. We employ 1,000 linear elements to 
resolve the small spatial scales (S — 0.001), and this results in a time scale of O(10 -6 ). 
Thus, we set the time step of the fine-scale computation at St = 10~ 6 . The number of 
fine-scale computations within each global time step is rif — 1000, and the time step 
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of the coarse integration is chosen as At c = n c St with n c — 49, 000. Thus, the global 
time step is At = (n/ + n c )5t = 0.05. The polynomial orders are set at K = 5 and 
L = 5 for the restrictions in the random space and the physical space, respectively. 
For this application, we achieve computational speed-up of ~ 50, compared to the 
full-scale MC simulation. The random solution profiles in Figure 14.81 show again 
good agreement between the full-scale MCFEM and the multiscale computation. We 
remark that the quantification of the computational speed-up is problem dependent. 
In the diffusion problems considered here, such speed-up is smaller at the beginning 
of the computation due to the fast evolution of the solution. However, the speed-up 
is significantly larger once the initial transient is passed. 





Fig. 4.8. Solution profile at T = 1. Left: mean solution, Right: standard deviation. 



5. Summary. In this paper we present an equation- free multiscale algorithm for 
integrating unsteady diffusion problems in a random medium with small-scale spatial 
structures, e.g., short correlation length. The method is based on the assumption 
that although the individual realizations of the random solutions are characterized 
by small scales, their ensemble averages are much smoother, characterized by larger 
scales. This motivates the use a set of coarse-grained observation variables, which are 
based on the pointwise polynomial approximations of the fine-scale random solutions. 
Such coarse-grained variables are approximated accurately on a coarse mesh, and 
integrated in time with large time steps. An equation- free approach is employed for 
the coarse integration, as the explicit knowledge of the governing equations of the 
coarse variables is unavailable in closed form. Details of the multiscale method is 
presented, and its accuracy and efficiency are documented by numerical examples. In 
particular, we demonstrate that the present constructions of the restriction and lifting 
operators allow us to successfully reconstruct representative fine-scale solutions based 
only on the knowledge of the coarse solutions. Future work will include a complete 
analysis of error estimates of the current method and applications to more complicated 
systems. 

It is interesting that, in this approach, one creates a "wrapper" around an ex- 
isting direct detailed solver, using it as a black box; traditional stochastic Galerkin 
algorithms would require the writing of new code to solve the coupled system of equa- 
tions in both physical and random space. Our approach can thus be considered as a 
"nonintrusive" one - the solution in both physical and random space is solved using 
an existing legacy code through a wrapper, and sidestepping the effort of new code 
development and validation. In this paper the only numerical task we demonstrated 
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in the equation-free context was temporal integration. Other tasks enabled through 
matrix-free iterative linear algebra (e.g. Newton-Krylov GMRES based fixed point 
solvers, Arnoldi-type eigensolvers) naturally fit in the equation-free framework; and 
many more can be performed on the explicitly unavailable coarse-grained equation: 
steady state and bifurcation computations, stability computations, equation-free op- 
timization and even dynamic renormalization. Finally, the smoothness in space of 
the coarse-grained variables can be exploited (via the so-called gap-tooth and patch 
dynamics equation free schemes) to perform the fully resolved fine-scale computations 
not only for short times, but also for only parts of the physical domain of interest 
[211 1361 IT^I . Finally, although in this paper we have designed direct numerical ex- 
periments to solve the coarse-grained equation which is hypothesized to exist and 
close, it is worth noting that it is also possible to design direct numerical experiments 
to test this hypothesis (see |25l 'l. These tasks, and the conditions under which they 
can be successfully performed in an equation-free framework and accelerate random 
computations, is the subject of ongoing research. 
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